
%% initialize
clear all;
close all;

rng(328120)

Nsim = 500;

loadstarting = 0;
loadestimates = 0;

%% read in data

% choices

rawdata = csvread('../../01_data/02_basedata/estimationdata/All_T_J/All_T_J_choices_buff_0_0_lf_tbar_Ndis.csv',1,0);

IJitidx         = rawdata(:,1);
IJjtidx         = rawdata(:,2);
IJapplied       = rawdata(:,3);
IJinterview     = rawdata(:,4);
IJchoiceset     = rawdata(:,5);
IJiidx          = rawdata(:,6);
IJjidx          = rawdata(:,7);
IJdistij        = rawdata(:,8);
IJteachidx      = rawdata(:,9);
IJschoolidx     = rawdata(:,10);
IJappyear       = rawdata(:,11);
IJcommute       = rawdata(:,15);
IJcommuteMiss   = rawdata(:,16);
IJhired         = rawdata(:,17);
IJofferdeclined = rawdata(:,18);
IJwithdrew      = rawdata(:,19);
IJpositive      = rawdata(:,20);
IJnegative      = rawdata(:,21);
IJVAma          = rawdata(:,22);
IJVAmean        = rawdata(:,23);
IJVAm1N          = rawdata(:,24);
IJVAm2N          = rawdata(:,25);
IJclassSize     = rawdata(:,26);
IJX             = rawdata(:,27:end-7);

IJVAmadev       = rawdata(:,27);
IJfracBlack     = rawdata(:,31);
IJfracHisp      = rawdata(:,32);
IJNDisadv       = rawdata(:,33);
IJmeanScore     = rawdata(:,34);

IJfracDisadv = IJNDisadv./IJclassSize;

IJX = [IJcommute,IJcommuteMiss,IJX,IJVAm1N,IJVAm2N];

IJcurrent       = rawdata(:,end-5);
IJcurrentFirst  = rawdata(:,end-4);
IJVAm1          = rawdata(:,end-3);
IJVAm2          = rawdata(:,end-2);


clearvars rawdata;

%% prepare data for estimation

IJchoiceset(IJcurrent==1) = 0;

d1 = max(IJitidx);
d2 = max(IJjtidx);

mint = min(IJappyear);
maxt = max(IJappyear);

IJcommutetrunc = min(35,IJcommute);


%% draw random effects

drawsi = normrnd(0,1,[d1,Nsim]);
drawsVAcomp = normrnd(0,1,[d1,Nsim]);
drawsBlack  = normrnd(0,1,[d1,Nsim]);
drawsHisp   = normrnd(0,1,[d1,Nsim]);
drawsDisadv = normrnd(0,1,[d1,Nsim]);
drawsScore  = normrnd(0,1,[d1,Nsim]);


Xdraws = zeros(size(IJitidx,1),6,Nsim);
Xdraws(:,1,:) = drawsi(IJitidx,:);
Xdraws(:,2,:) = drawsi(IJitidx,:).*(IJVAmadev*ones(1,Nsim));
Xdraws(:,3,:) = drawsi(IJitidx,:).*(IJNDisadv*ones(1,Nsim));
Xdraws(:,4,:) = drawsVAcomp(IJitidx,:).*(IJVAmadev*ones(1,Nsim));
Xdraws(:,5,:) = drawsVAcomp(IJitidx,:).*(IJNDisadv*ones(1,Nsim));
Xdraws(:,6,:) = drawsDisadv(IJitidx,:).*(IJNDisadv*ones(1,Nsim));


%% estimate mixed logit


y_logistic = IJapplied(IJchoiceset==1);

   
Xmat = [ones(length(IJapplied),1),IJX];

f =  @(x)func_mixed_logit_current(x,IJapplied(IJchoiceset==1),Xmat(IJchoiceset==1,:),...
    Xdraws(IJchoiceset==1,:,:),IJitidx(IJchoiceset==1),Xmat(IJcurrentFirst==1,:),...
    Xdraws(IJcurrentFirst==1,:,:),IJitidx(IJcurrentFirst==1),0);


    betastart = zeros(size(Xmat,2),1);
    sigmastart  = ones(size(Xdraws,2),1); sigmastart(3) = 1/60; sigmastart(5) = 1/60; sigmastart(6) = 1/60;
 
    thetastart = [betastart;sigmastart];





optionsgrad = optimset('MaxFunEvals',10000000,'MaxIter',10000000,'TolX',1E-5,...
    'Display','Iter','TolFun',1E-7,'Algorithm','trust-region','GradObj','on');


    tic
thetahat = fminunc(f,thetastart,optionsgrad);
toc
f =  @(x)func_mixed_logit_current(x,IJapplied(IJchoiceset==1),Xmat(IJchoiceset==1,:),...
    Xdraws(IJchoiceset==1,:,:),IJitidx(IJchoiceset==1),Xmat(IJcurrentFirst==1,:),...
    Xdraws(IJcurrentFirst==1,:,:),IJitidx(IJcurrentFirst==1),1);
[obj,grad,se,p1mat] =  f(thetahat);

 save -v7.3 '../../01_data/05_temp/mixed_logit_rc6_totVA_DISAD_corrRC'
 
estimates_se = [thetahat,se];
writematrix(estimates_se,'../../03_output/tables/teacher_pref_estimates_rc6_totVA_DISAD_corrRC.csv')

